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Abstract 

We analyze the random fluctuations of several multi-scale algorithms such as the 
multi-scale finite element method (MsFEM) and the finite element heterogeneous multi- 
scale method (HMM), that have been developed to solve partial differential equations 
with highly heterogeneous coefficients. Such multi-scale algorithms are often shown to 
correctly capture the homogenization limit when the highly oscillatory random medium 
is stationary and ergodic. This paper is concerned with the random fluctuations of 
the solution about the deterministic homogenization limit. We consider the simplified 
setting of the one dimensional elliptic equation, where the theory of random fluctuations 
is well understood. We develop a fluctuation theory for the multi-scale algorithms in 
the presence of random environments with short-range and long-range correlations. 

For a given mesh size h, we show that the fluctuations converge in distribution in 
the space of continuous paths to Gaussian processes as the correlation length e —> 0. 
We next derive the limit of such Gaussian processes as h — > and compare this limit 
with the distribution of the random fluctuations of the continuous model. When such 
limits agree, we conclude that the multi-scale algorithm captures the random fluctua- 
tions accurately and passes the corrector test. This property serves as an interesting 
Benchmark to assess the behavior of the multi-scale algorithm in practical situations 
where the assumptions necessary for the theory of homogenization are not met. 

What we find is that the computationally more expensive methods MsFEM, and 
HMM with a choice of parameter 5 = h, correctly capture the random fluctuations 
both for short-range and long-range oscillations in the medium. The less expensive 
method HMM with 5 < h correctly captures the fluctuations for long-range oscillations 
and strongly amplifies their size in media with short-range oscillations. We present a 
modified scheme with an intermediate computational cost that captures the random 
fluctuations in all cases. 

Keywords: Equations with random coefficients, multi-scale finite element method, 
heterogeneous multi-scale method, corrector test, correlation ranges. 

AMS subject classification: 65C99, 65Y20, 34E13 



1 Introduction 

Differential equations with highly oscillatory coefficients arise naturally in many areas of 
applied sciences and engineering from the analysis of composite materials to the modeling of 
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geological basins. Often, it is impossible to solve the microscopic equations exactly, in which 
case we can either solve macroscopic models when they exist or devise multi-scale algorithms 
that aim to capture as much of the microscopic scale as possible. When the coefficients are 
periodic functions or stationary and ergodic processes, the solution to the heterogeneous 
equation is often well approximated by the deterministic solution to a homogenized equation. 
This is the well known homogenization theory. In many applications, such as parameter 
estimation and uncertainty quantification, estimating the random fluctuations (finding the 
random corrector) in the solution is as important as finding its homogenized limit. 

Finding the homogenized coefficients may be a daunting computational task and the 
assumptions necessary to the applicability of homogenization theory may not be met [18, 
20, 22]. Several numerical methodologies have been developed to find accurate approxi- 
mations of the solution without solving all the details of the micro-structure [1, 10, 12]. 
Examples include the multi-scale finite element method (MsFEM) and the finite element 
heterogeneous multi-scale method (HMM). Such schemes are shown to perform well in the 
homogenization regime, in the sense that they approximate the solution to the homogenized 
equation without explicitly calculating any macroscopic, effective medium, coefficient. Ho- 
mogenization theory thus serves as a Benchmark that ensures that the multi-scale scheme 
performs well in controlled environments, with the hope that it will still perform well in 
non-controlled environments, for instance when ergodicity and stationarity assumptions are 
not valid. 

This paper aims to present another Benchmark for such multi-scale numerical schemes 
that addresses the limiting stochasticity of the solutions. We calculate the limiting (prob- 
ability) distribution of the random corrector given by the multi-scale algorithm when the 
correlation length of the medium tends to at a fixed value of the discretization size h. We 
then compare this distribution to the distribution of the corrector of the continuous equa- 
tion. When these distributions are close, in the sense that the h— dependent distribution 
converges to the continuous distribution as h — > 0, we deduce that the multi-scale algorithm 
asymptotically correctly captures the randomness in the solution and passes the random 
corrector test. 

The above proposal requires a controlled environment in which the theory of correctors 
is available. There are very few equations for which this is the case [2, 3, 14]. In this paper, 
we initiate such an analysis in the simple case of the following one-dimensional second-order 
elliptic equation: 



Here, the diffusion coefficient a(|,u>) is modeled as a random process, where u denotes the 
realization in an abstract probability space (f2, J-, P) in which the random process and all 
limits considered here are constructed. The correlation length e is much smaller than L = 1, 
the length of the domain, which makes the random coefficient highly oscillatory. 

We chose this equation for two reasons. First, many multi-scale numerical schemes 
have been developed to solve its generalization in higher dimensions. Second, both the 
homogenization and the corrector theory for this elliptic equation are well known. When 
a(x,co) is stationary, ergodic, and uniformly elliptic, i.e., A < a(x,co) < A for almost 
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every x and u, then the solution u e converges to the following homogenized equation with 
deterministic and constant coefficient (this generalizes to higher dimensions as well [20, 22]): 

uo(0) =«o(l) = 0. 

The coefficient a* is the harmonic mean of a(x, lj). Furthermore, if the deviation of l/a(x, oj) 
from its mean 1/a* is strongly mixing with integrable mixing coefficient, a notion that will 
be defined in the next section and which in particular implies that the random coefficient 
has short range correlation (SRC), then the corrector theory in [6] shows that 

u e ~ u I \ distribution „t\ /-. „x 

(x) >U(x;W), (1.3) 



where U{x\ W) is a Gaussian random process that may be conveniently written as a Wiener 
integral; see (2.13) below. The convergence above is in the sense of distribution in the space 
of continuous paths, which will be denoted by C([0, 1]). When the deviation 1/a — 1/a* 
does not de-correlate fast enough, the normalization factor yfe is no longer correct. In 
fact, for a large class of random processes with long range correlation (LRC) defined more 
precisely later, it is proved in [3] that a similar convergence result holds with the following 
modifications: The normalization factor should be replaced by e'a, where < a < 1 
describes the rate of de-correlation. The limiting process is then Ujj(x;W h ), which is 
defined in (2.14) and has the form of a stochastic integral with respect to fractional Brownian 
motion. 

Let us assume that we want to solve (1.1) numerically. We denote by h the discretization 
size and by u^(x) the solution of the scheme. The standard finite element solution of the 
deterministic limit (1.2) is denoted by Uq(x). Our goal is to characterize the limiting 
distribution of random process after proper rescaling by e 2 . Here and below, 

we use c Ad to denote min{c, d}. Obviously, this normalization parameter is chosen to be 
consistent with the aforementioned corrector theories. 

We say that a numerical procedure is consistent with the corrector theory and that it 
passes the corrector test when the following diagram commutes: 

u e~ u 0f.. \ h->0. Ue-UQ, 



— (a?,w) — — ■> ^—{X,U)) 



aAl v • , ... 
£2 W £2 



(ii) (iii) 



£^0 



;i-4) 



U^ M (x;W aM ) U aM (x;W c 

(tv) 



Here, W 1 = W is the standard Brownian motion whereas W a = W H is the fractional 
Brownian motion with Hurst index H = 1 — ^. Similarly, U a /\i is defined so that IA\ = 
U(x;W) and U a = Uh{x;W h ). We use this notation to include both random processes 
with SRC and those with LRC. In the above diagram, there are four convergence paths 
to be understood in the sense of distribution in C([0, 1]). The convergence path (iii) is 
the corrector theories in [6, 3]. In (i), h is sent to zero while e is fixed. To check (i) is a 
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numerical analysis question without multi-scale issues since the e-scale details are resolved. 
Convergence in (i) can be obtained path-wise and not only in distribution (path (iv) may 
also be considered path- wise). The main new mathematical difficulties we address in this 
paper, therefore lie in analyzing the paths (ii) and (iv). 

Our main results are stated for the multi-scale finite element method (MsFEM) [16, 17, 
15], and the finite element heterogeneous multi-scale method (HMM) [9, 10]. The main 
idea of MsFEM is to use multi-scale basis functions constructed from the local solutions of 
the elliptic operator in (1.1) on domains of size h <C 1. The main idea of HMM is to con- 
struct these basis functions on possibly smaller patches. HMM thus involves heterogeneous 
computations on domains of size 5 < h. 

The analysis of MsFEM in random settings is done in [7, 11]. HMM in random media 
was considered in [10]. These works show that both MsFEM and HMM pass the test of 
homogenization theory in the sense that as e — > 0, the limiting solution is a consistent 
discretization (i.e., with error converging to as h — > 0) of the homogenized limit. This 
paper analyzes the corrector theory of both multi-scale methods. 

For MsFEM and HMM with the parameter choice 5 = h (these two methods are then 
the same in dimension n = 1), the above diagram commutes. Path (i) holds as a result of 
standard finite element analysis. Path (ii) holds due to the self- averaging effect of integrals 
with oscillatory integrand a~ 1 (x/e)—a*~ 1 , andU^ Al is explicitly characterized as a Gaussian 
process. The expression for this process is complicated but converges to the right corrector 
as h — > 0. 

For HMM with 5 < h, path (i) does not hold because there is no homogenization effect 
when h <C £■ Passing the corrector test then means that (ii) followed by (iv) yields the 
same result as (iii). Path (ii) is obtained with U^ A1 a Gaussian process. HMM with 6 < h 
then performs differently for SRC and LRC. In the first case, the limit process as h — > in 
(iv) is {h/5)^TA: the limit is an amplified version of the theoretical corrector when 5 < h. 
In that sense, HMM with 5 < h does not pass the corrector test. In the LRC case, however, 
and somewhat surprisingly at first, HMM passes the corrector test for all < 5 < h. 
As a consequence, HMM with 5 < h finds the corrector up to a rescaling coefficient that 
depends on the structure of the random medium and whose estimation may prove difficult 
for environments in which the corrector theory is not valid or not known. 

We consider a modification of the MSFEM and HMM schemes that is computationally 
more expensive than HMM with 5 <C h but passes the corrector test independent of the 
underlying random structure. The proposed corrector test and our analysis, in spite of 
the limited scope of the one-dimensional equation, thus provide some guidance on the 
accuracy one may expect from a multi-scale algorithm based on a given computational cost 
and a given level of prior information about the heterogeneous medium (such as, e.g., its 
correlation properties) one is willing to make. 

The rest of the paper is organized as follows. In the next section, we present our models 
for the random medium and formulate the main results of the paper. The derivations rely 
on establishing explicit expressions for — Uq . In section 3, we derive such expressions for 
a general numerical scheme satisfying specific assumptions on the structure of the associated 
stiffness matrix. We show how to prove (ii) using these formulas. In section 4, we show how 
to prove (iv) in the general setting. Both convergences in (ii) and (iv) can be viewed as weak 
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convergence of measures in C, which can be shown as usual by obtaining the convergence 
of finite dimensional distributions and establishing tightness. In section 5, we apply the 
general framework to MsFEM and prove that diagram (1.4) commutes in this case. In 
section 6, we consider HMM in both media with SRC or LRC. We emphasize where the 
amplification effect comes from in the case of SRC and why this effect disappears in the case 
of long-range media. In section 7, we discuss methods to eliminate the amplification effect 
of HMM. Some conclusions are offered in section 8 while technical lemmas on fourth-order 
moments of random processes are postponed to the appendix. 

2 Main results on the corrector test 

In this section, we introduce our hypotheses on the random media, describe the multi-scale 
algorithms MsFEM and HMM and formulate our main convergence results. 

To simplify notation, we drop the dependency in oj when this does not cause confusion. 
We define a e (x) = a(x/e). For a function g in LP(D), we denote its norm by ||g||p,D- 
When D is the unit interval we drop the symbol D. The natural space for (1.1) and 
(1.2) is the Hilbert space H^. By the Poincare inequality, the semi-norm of Hq defined by 
\ u \h 1 ,d = H^A^Ih.-Dj is equivalent to the standard norm. We use the notation C to denote 
constants that may vary from line to line. When C depends only on the elliptic constants 
(A, A), we refer to it as a universal constant. Finally, the Einstein summation convention 
is frequently used: two repeated indices, such as in CjoP, are summed over their (natural) 
domain of definition. 

2.1 Random media models 

We model a(x,co) as a random process on some probability space (fl, P), with parameter 
iGR. The coefficient a(|,w) in (1.1) is obtained by rescaling the spatial variable. 

A process a(x) is called stationary if the joint distribution of any finite collection of 
random variables {a(xj, w)}i<j< n . is translation invariant. Let E denote the mathematical 
expectation with respect to P. We define 



Since a is stationary, a* is a constant and q(x,u) is a stationary process with mean zero. 
The (auto-)correlation function R(x) of q is given by: 



We find that R is symmetric, i.e., R(x) = R(—x) and is a function of positive type in the 
sense of [25, §XI]. By Bochner's theorem, the Fourier transform of R is non-negative. Let 
us define 




a* 



1 



(2.1) 



R(x) := E{q(y)q(y + x)} = K{q(0)q(x)}. 



(2.2) 




(2.3) 
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Then a 2 is in [0, oo], and we assume here that a > 0. We say that q(x,co) has short range 
correlation (SRC) if R is integrable. In this case, a is a positive finite number. We say q 
has long range correlation (LRC) if the correlation function is not integrable. 

The decay of R(x) at infinity describes the two-point decorrelation rate of the process 
q. We also need the notion of mixing. Consider a Borel set A C K. Let J- a denote the 
Borel sub-cr-algebra of J- induced by {q(x);x G ^4}, and L A the space of J- a— measurable 
and square integrable random variables. A random process a(x,u)) is said to be p-mixing 
if there exists some non-negative function p(r) : M+ —> R+ with p(r) — > at infinity, such 
that for any Borel sets A and B the following holds. 

sup EM-ssffijU (2 . 4) 

teL% v eL% Var{^}2Var{r/}2 

Here d(^4, B) is the distance between the two sets. The function p(r) is called the /j-mixing 
coefficient, and it reflects, roughly speaking, how fast random processes restricted on A and 
B become independent. More details on the notion of mixing can be found in [8]. 

Our first set of assumptions on a(x) and q(x) are the following. We note that the third 
assumption below implies the second one. 

(51) The random process a{x) is stationary and uniformly elliptic with constants (A, A). 
Clearly, a* is uniformly elliptic with the same constants. 

(52) The random process q{x) has SRC and we assume a > in (2.3). 

(53) The mixing coefficient p(r) of q(x) is decreasing in r, and p5( r ) g L 1 (M + ). 

The above assumptions are quite general. In particular, (S3) implies ergodicity of q(x), 
and (SI) plus ergodicity is the standard assumption for homogenization theory; (S3) is the 
standard assumption to obtain central limit theorem of oscillatory integrals with integrand 
q £ (x) as in (3.29). 

For random media with LRC, there is no general central limit theorem. A special family 
of random media investigated in [3] is based on the following assumptions: 



(LI) The process q(x) is constructed as 

q(x,U)) = $(g(x,cj)), (2.5) 

where g x is a stationary Gaussian random process with mean zero and variance one. 
Further, the correlation function R g of g x has the following asymptotic behavior: 

Rg(r) ~ K g r- a , (2.6) 

where n g > is a constant and a G (0, 1). 

(L2) The function <3?(x) satisfies |$| < go f° r some constant qo, so that the process a(x,e), 
constructed by the relation (2.1) for some positive constant a* , satisfies uniform ellip- 
ticity with constants (A, A). 



G 



(L3) The function $ integrates to zero against the standard Gaussian measure: 



$(s)e~irds = 0. (2.7) 

The process q(x) is stationary and mean-zero. More importantly, its correlation function 
R{x) has a similar asymptotic behavior to that in (2.6) with K g replaced by k, where 

1 r s 2 \ 2 

s<S>(s)e~^ds . (2.8) 



The integral above is assumed to be non-zero. Consequently, R(x) is no longer integrable 
and q(x) has LRC. We note that when a > 1, the process constructed above has SRC and 
provides an example satisfying (S2). 

2.2 The multi-scale algorithms MsFEM and HMM 

We briefly introduce MsFEM and HMM and leave some details to later sections. Assume 
/£l 2 C H -1 . The weak solution to (1.1) is the unique function u £ S -£^(0, 1) such that 

A £ (u,v) = F(v), Vu 6^(0,1). (2.9) 

The associated bilinear and linear forms are defined as 

A £ (u,v) : = J\ £ (x)^ ~dx, F(v) := J fv dx. (2.10) 

Existence and uniqueness of u e are guaranteed by the uniform ellipticity of a e (x). MsFEM 
and the HMM are finite element methods and so the ideas are to approximate Hq by finite 
dimensional space and, when necessary, to approximate A e by an auxiliary bilinear form. 

We partition the unit interval into N small intervals of size h = 1/N. We use this 
partition for MsFEM, HMM, and the standard FEM. The FEM basis functions are piece- 
wise linear "hat functions" , each of them peaking at one nodal point and vanishing at all 
other nodal points. Denote these hat functions by {<Pq} and denote the subspace of Hq they 
span by V h . The standard FEM approximates Hq by V h . The idea of MsFEM is to replace 
the hat basis functions by multi-scale basis functions {4>l}- For instance, (fP £ is constructed 
by solving the elliptic equation (1.1) locally on the support of cfy; see (5.1) below. Denote 
the span of these multi-scale basis functions by V £ . Then, the MsFEM solution, denoted 
by u £ , is the unique function in V £ so that (2.9) holds for all v S V £ . The well-posedness 
of this weak formulation is a consequence of the uniform ellipticity of a £ . 

HMM aims to reduce the cost of MsFEM by computing heterogeneous solutions on 
smaller domains. Solutions and test functions are still sought in V h . The bilinear form A £ 
is then approximated by 




A h /(u, v) := h / . a ^^T " dx~^ dX K 
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where it is a small sub-interval of size 5 < h of the A;th interval in the above partition. The 
operator is defined below in (6.1). This linear operator optimally uses the micro-scale 
calculation on it in the same way that (f>i does in MsFEM. HMM solution is the unique 
function u £ ' S in Vq so that (2.9) holds when Hq is replaced by V h and A £ is replaced by 
A £ ' S . The well-posedness of this formulation is less immediate and will be proved below. 

Throughout this paper, for simplicity in presentation, we assume that the micro-scale 
solvers in MsFEM and HMM, i.e., equations (5.1) and (6.1), are exact. 



2.3 Main theorems 



We state our main convergence theorems in the setting of MsFEM and HMM although they 
hold for more general schemes. 

The first theorem analyzes MsFEM in the setting of SRC. 

Theorem 2.1. Let u £ and uq be solutions to (1.1) and (1.2), respectively. Let u £ be the 
solution to (1.1) obtained by MsFEM and letu^ be the standard finite element approximation 
of uq. Then we have: 

(i) Suppose that a(x) satisfies (SI) and f is continuous on [0, 1]. Then 

\u h £ - u £ \ m < An/112, \\ u h _ U£ || 2 < A^||/|| 2 . (2.11) 



Assume further that q(x) satisfies (S2). Then 



sup 

x€[0,l] 



<C^||i?|| 1>M (l + ||/|| 2 ), (2.12) 



where C is a universal constant and R is the correlation function of q as defined in (2.2). 
(ii) Now assume further that q(x) satisfies (S3). Then, 

distribution / 



>a [ L h (x,t)dW t =:U h {x;W). (2.13) 
Jo 



The constant a is defined in (2.3) and W is the standard Wiener process. The function 
L h (x,t) is explicitly given in (3.20). The convergence is in distribution in the space C. 

(iii) Now let h goes to zero, we have 

U h {x; W) dlstnbutlon > u(x; W):=a / L{x,t)dW t . (2.14) 
h->0 Jo 

The Gaussian process U{x\ W) was characterized in [6]. The kernel L(x,t) is defined as 

L(M) = l[o,*](*)(7 F( 8 )ds-F(t)j+x(F(t)-J Q F(s)ds^jl m (t). (2.15) 
Here and below, 1 is the indicator function and F(t) = f Q l f(s)ds. 



S 



Equivalently, the theorem says the diagram in (1.4) commutes when q has SRC and that 
MsFEM passes the corrector test in this setting. 

To prove (iv) of the diagram, we recast L(x, t) as 

L(x, t) = a*-^G (x, t) ■ a*foMt)- ( 2 - 16 ) 
Here Go is the Green's function of (1.2). It has the following expression: 

G (x,y) = l ^_ (2.17) 
[a (l-x)y, x>y. 

In particular, Go is Lipschitz continuous in each variable while the other is kept fixed. 

The next theorem accounts for MsFEM in media with LRC. We recall that the random 
process q(x) below is constructed as a function of a Gaussian process. 

Theorem 2.2. Let u £ , Uq, u £ and Uq be defined as in the previous theorem. Let q(x,u) 
and a(x,u) be constructed as in (L1)-(L3). Then we have 

(i) 

E(u h £ (x)-u h (x)) 



sup 

xe[o,i] 



2 



<cUiY, (2.18) 



h \h 



for some constant C depending on (A, A), k, a, and f. 
(ii) As e goes to zero while h is fixed, we have 

^o(^) distribution 7 ./j 



£2 £->0 



>U^(x;W H ):=a H [ L h (x,t)dW t H . (2.19) 



Here H = 1 — ^ , and W^ 1 is the standard fractional Brownian motion with Hurst index H. 
The constant o~h is defined as n/H(2H — 1). The function L h (x,t) is defined as in the 
previous theorem. 

(iii) As h goes to zero, we have 

U^{x;W H ) distributi011 ) U H (x;W H ) := a H C L(x,t)dW t H . (2.20) 

As before, this theorem says the diagram in (1.4) commutes in the current case. In 
particular, a < 1, and the scaling is £2 . Thus MsFEM passes the corrector test for both SRC 
and LRC. The stochastic integrals in (2.19) and (2.20) have fractional Brownian motions 
as integrators. We give a short review of such integrals in the appendix. A good reference 
is [23]. 

The next theorem addresses the convergence properties of HMM. 

Theorem 2.3. Let u e and uq be the solutions to (1.1) and (1.2), respectively. Let u £ ' S be 
the HMM solution and Uq the standard finite element approximation ofuQ. 
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(i) Suppose that the random processes a(x) and q(x) satisfy (S1)-(S3). Then 

U h £ '\x)-U h (x) distribution ) ^ , _ distribution ) Jh 
yje e-s-0 /i-vO V o 

Here, U h ' S (x;W) is as in (2.13) with L h replaced by L h,s (x,t), which is defined in (6.7) 
below. The process U(x;W) is as in (2.14). 

(ii) Suppose instead that the random processes a(x) and q(x) satisfy (L1)-(L3). Then 

U*'\x)-V%(x) distribution U H 5, X . wH) distribution U H (x] W H ). (2.22) 

Here, U^ 5 (x;W H ) is as in (2.19) with L h replaced by L h > 5 , andU H {x;W H ) is as in (2.20). 

HMM is computationally less expensive than MsFEM when 5 is much smaller than 
h. However, the theorem implies that this advantage comes at a price: when the random 
process q(x) has SRC, HMM with 5 < h amplifies the variance of the corrector. We will 
discuss methods to eliminate this effect in section 7. In the case of LRC, however, HMM 
does pass the corrector test even when 6 <C h. 

Intuitively, averaging occurs at the small scale 5 <C h for SRC. Since HMM performs 
calculations on a small fraction of each interval h, each integral needs to be rescaled by h/5 
to capture the correct mean, which over- amplifies the size of the fluctuations. In the case 
of long-range correlations, the self-similar structure of the limiting process shows that the 
convergence to the Gaussian process occurs simultaneously at all scales (larger than e) and 
hence at the macroscopic scale. HMM may then be seen as a collocation method (with grid 
size h), which does capture the main features of the random integrals. 

The amplification of the random fluctuations might be rescaled to provide the correct 
answer. The main difficulty is that the rescaling factor depends on the structure of the 
random medium and thus requires prior information or additional estimations about the 
medium. For general random media with no clear scale separation or no stationarity as- 
sumptions, the definition of such a rescaling coefficient might be difficult. In section 7, we 
present a hybrid method between HMM and MsFEM that is computationally less expensive 
than the MsFEM presented above while still passing the corrector test. 



3 Expression for the corrector and convergence as e — > 

The starting point to prove the main theorems is to derive a formula for the corrector u^—Uq. 
This is done by inverting the discrete systems yielded from the multi-scale schemes. The 
goal is to write the corrector as a stochastic integral linear in the random coefficient, plus 
error terms that are negligible in the limit; see Proposition 3.2 below. 



3.1 General finite element based multi-scale schemes 

Almost all finite element based multi-scale schemes for (1.1) have the same main premise: 
in the weak formulation (2.9), we approximate Hq by a finite dimensional space and if 
necessary, also approximate the bilinear form. 
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To describe the choices of the finite spaces, we choose a uniform partition of the unit 
interval into N sub-intervals with size h = 1/N. Let x k denote the kih. grid point, with 
xq = and xjy = 1, and I k the interval (x k _i,x k ). To simplify notation in the general 
setting, we still denote by V £ the finite space and by {(ftl}^^ 1 the basis functions. They do 
not necessarily coincide with those in MsFEM. 

In a general scheme, the bilinear form in (2.9) might be modified. Nevertheless, to 
simplify notations, we still denote it as A £ . The solution obtained from the scheme is then 
u £ G V £ such that 

A £ {u h £ ,v) = F{v), Vv^V £ h . (3.1) 

Since V £ is finite dimensional, the above condition amounts to a linear system, which is 
obtained by putting u £ = Uj(j){, and by requiring the above equation to hold for all basis 
functions. The linear system is: 

A h v e = F e (3.2) 

Here, the vector U £ is a vector in R^ -1 , and it has entries t/J. The load vector F e is also 
in IR^" 1 and has entries F((f) £ ). The stiffness matrix A £ is an N — 1 by N — 1 matrix, and 
its entries are A £ ((f) l £ ,(jP £ ). Our main assumptions on the basis functions and the stiffness 
matrix are the following. 

(Nl) For any 1 < j < N — 1, the basis function (pi is supported on Ij U Ij+i, and it takes 
the value 5j at nodal points {x{\. Here 5j is the Kronecker symbol. 

(N2) The matrix A £ is symmetric and tri-diagonal. In addition, we assume that there exists 
a vector b £ G M. N with entries {b £ }^ =l , so that A £ii+1 = — b £ +1 for any i = 1, • • • , N — 2 
and 

A h £ii = -(A h £ii _ 1 + A h £ii+1 ), i = l,---,N-l. (3.3) 

In other words, the ith. diagonal entry of A £ is the negative sum of its neighbors in 
each row. Here, ^4^01 ^rid A £ j^_ -^jy are not matrix elements and are set to be b £ and 
b^ , respectively. 

(N3) On each interval for j = 1, • • • , N, the only two basis functions that are nonzero 
are cpi and (fP e 1 , and they sum to one, i.e., <fp e + 1 = 1. Equivalently, we have 

4(x) = + - & + \x)], (3.4) 

for some functions {4> £ {x)}^ =l , each of them defined only on Ij with boundary value 
at the left end point and 1 at the right. 

As we shall see for MsFEM, (N3) implies (N2) when the bilinear form is symmetric. 
The special tri-diagonal structure of A £ in (N2) simplifies the calculation of its action on a 
vector. Let U be any vector in R N_1 , we have 

(A h £ U) i = -D + (b i £ D-U i ), i = !,-■■ ,N-1. (3.5) 

Here, the symbol D~ denotes the backward difference operator, which is defined, together 
with the forward difference operator D + , as 

(D-U) k = U k -U k . u (D + U) k = U k+l -U k . (3.6) 
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The equality (3.5) is easy to check, and to make sense of the case when i equals 1 or N, 
we need to extend the definition of U by setting Uq and Un to zero. This formula has been 
used, for example, in [17]. It is a very useful tool in the forthcoming computations. 



3.2 The expression of the corrector 

Now we derive an expression of the corrector, i.e., the difference between the solution 
to (1.1) obtained from the above scheme, and Uq, the standard FEM solution to (1.2). 

The function Uq(x) is obtained from a weak formulation similar to (2.9) with a £ replaced 
by a*, and Hq replaced by V h , the space spanned by hat functions {<ft }. Clearly, these 
basis functions satisfy (Nl) and (N3). Let Aq denote the associated stiffness matrix; then 
one can verify that it satisfies (N2). In fact, the vector b is given by 6q = a*/h. Now Uq(x) 
is simply U^cfy, where U° solves 

A%U° = F°. (3.7) 
Subtracting this equation from (3.2), we obtain: 

A%(U £ - U°) = {F £ - F°) - {A h £ - A^)U £ . 

Let Gq denote the inverse of the matrix Aq. We have 

U £ -U° = G^{F £ - F°) - G^{A h e - A%)U £ . 

Since both A^. and A^ satisfy (N2), the difference A^. — A^ acts on vectors in the same 
manner as in (3.5). Since both {0e} and {4>q} satisfy (N3), we verify that 

{F e _ F 0) . = _ D+[ pe _ pO^ pe . = f f{t) ^ [t)dL 

Jlj 

Using these difference forms, we have 



N-l 

U J- U J=-H (d + (F £ - F°) m - D + ((bf - b%)D-U £ m 

m=l 
N 

= £ D~G h 0jk ((F £ - F\ - (b k £ - b k )D~U £ k 
k=l 



(3.8) 



The second equality is obtained from summation by parts. Note that we have extended the 
definitions of U £ and U° so that they equal zero when the index is or N. Similarly, (Gq)jo 
and {Gq)jn are zero as well. 

The vector U £ — U° is the corrector evaluated at the nodal points. We have the following 
control of its I 2 norm under some assumptions on the statistics of {b k } and {4>i}. 

Proposition 3.1. Let U £ and U° be as above. Let the basis functions {4>i} an d the stiffness 
matrix A^ satisfy (N1)-(N3). Suppose also that 

sup \D-U £ k \<C\\f\\ 2 h l 2, (3.9) 

Kk<N 
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for some universal constant C . 

(i) Suppose further that for any k 



,N, and any x E Ik, we have 



E $(3) ) <C-\\R\\ K 



V(b k £ -b k Q ) ^C^WRh 



and 

for some universal constant C . Then we have 

K\\U £ 



U°\\l <Cp||i?|| 1)R (l + 



(3.10) 
(3.11) 

(3.12) 



for some universal C . 

(ii) Suppose instead that the right hand side of (3.10) is C (§)", and the right hand side 
of (3.11) is Cjp (f)°- Then the estimate in (3.12) should be changed to C-r? (§)". 

The assumption (3.9) essentially says that should have a Holder regularity. Suppose 
the weak formulation associated to the multiscale scheme admits a unique solution such 
that \\u^ \\ H i < C||/||2. Then by Morrey's inequality [13, p. 266], G C°'5 in one dimension. 
Consequently, (3.9) holds. 

For MsFEM, we have a better estimate: \D~U £ \ < Ch due to a super-convergence 
result; see (5.3). Therefore, the estimate in (3.12) can be improved to be C-rj in case (i) 
and Cj^ (§)° in case (ii). 

Item (i) of this proposition is useful when the random medium a(x), or equivalently 
q(x), has SRC, while item (ii) is useful in the case of LRC. The constant C in the second 
item depends on (A, A), / and R g , but not on h. 



Proof. To prove (i), we use a super-convergence result, which we prove in section 5.1, to 
get \D~GQj k \ < Ch. Using this estimate together with (3.9) and (3.8), we have 

N _ N 



E|[/J 



U 



0|2 



<ChJ2 E \Fk-Fk +C7i 2 ^E^-6 l 



k=l 



k=l 



For the second term, we use (3.11) and obtain 

N N 



^Ek fc - b k " <^C^\\R\\ hR = C^\\R\\ X> 



k=l 



k=l 



(3.13) 



For the other term, an application of Cauchy-Schwarz to the definition of F £ yields 



Using (3.10), we have 



El 



h k \\ 2 



n 



E 



< 



T/c 1 1 2 



2 (x)dx <C^-h\\R\\ K 



Ce\\R\\i- 



(3.14) 
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Therefore, we have 

E \Uj - U°\ 2 < (chJ2\\fh,i k e + Ch 2 ^j \\R\\ 1M < C^||i?|| 1)K (l + ||/|| 2 ). 

Note that this estimate is uniform in j. Sum over j to complete the proof of (i). 

Proof of item (ii) follows in exactly the same way, using the corresponding estimates. □ 

Now, the corrector in this general multi-scale numerical scheme is: 

u h e {x) - u h (x) = U?4{x) - U?4(x) 

= {W - U%4(x) + U?(4 - 4)(x) + - U%(4 - 4)(x). 

We call the three terms on the right hand side Ki(x), i = 1,2,3. Now K\(x) is the 
piecewise interpolation of the corrector evaluated at the nodal points; ^(x) is the corrector 
due to different choices of basis functions; and K%(x) is much smaller due to the previous 
proposition and (3.10). Our analysis shows that K\{x) and K2(x) contribute to the limit 
when e — > while h is fixed, but only a part of K±(x) contributes to the limit when h — > 0. 

Due to self- averaging effect which is made precise in Lemma A.l, integrals of q e (x) 
are small. Therefore, our goal is to decompose the above expression into two terms: a 
leading term which is an oscillatory integral against q £ , and a remainder term which contains 
multiple oscillatory integrals. 

Proposition 3.2. Assume that u £ is the solution to (1.1) obtained from a multi-scale 
scheme, which satisfies (N1)-(N3) and has basis functions {4}; an d that Uq is the solution 
of (1.2) obtained by the standard FEM with hat basis functions {4q\- Let b e and bo denote 
the vectors in (N2) of these methods. Suppose that (3.9) holds and that for any k = 1, • • • ,N, 
we have 



l k £ {t) - 4(t) = [1 + f*]£ ( f Qe{s)ds - £ q £ (s)ds 



h k_ h k = [1 + hk] {_"_ j 



2 rx k 



(3.16) 



for some random variables fn- and r^. 

(i) Assume that q{x) has SRC, i.e., satisfies (S2), and that 

sup max{E|fi fc | 2 ,E|f 2fc | 2 } < \\R\\i,u, (3.17) 

l<fc<7V h 

for some universal constant C. Then, the corrector can be written as 

u h £ {x) - uft(x) = f L h (x,t)q £ (t)dt + r^(x). (3.18) 



o 



Furthermore, the remainder r^(x) satisfies 



sup E\r^(x)\ <C-^||fl||i i R(l + ||/|| 2 ), (3.19) 
xe[o,i] 
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for some universal constant C. The function L h (x,t) is the sum of L\ and L\ defined by: 



N 



L i(x,t) = 2^1j h (t) +/ f(s)ds- f(s)(j) {s)ds 

k=l \ Jt Jxk-i 

^(M) = U Oj(x) (ifaw-trtit) jp- MxjW-i^jix)}^ 

(3.20) 

Given x, the index j(x) is the unique one so that Xjr x \i < x < Xj( x y The function 
GQ(x,Xk) is defined as 

N-l 

G h (x,x k ) = Y,G h 0j Mx)- (3-21) 
i=i 

Gq is the interpolation in Vq using the discrete Green's function of standard FEM. 

(ii) Assume that q(x) has long range correlation, i.e., (L1)-(L3) are satisfied, and that 
the estimate in (3.17) is C (§)"• Then the same decomposition holds, the expression of 
L h (x,t) remains the same, but the estimate in (3.19) should be replaced by (§)"• 

Due to the super-convergent result in section 5.1, the function Gq(x, x k ) above is exactly 
the Green's function evaluated at (x,x k ). This can be seen from the facts that they agree 
at nodal points and both are piece-wise linear and continuous. 

Proof. We only present the proof of item (i). Item (ii) follows in exactly the same way. 
We point out that the assumption (3.16) and the estimates (3.17) imply (3.10) and (3.11) 
thanks to Lemma A.l. 

The idea is to extract the terms in the expression (3.15) that are linear in q e . For K\(x), 
we use (3.8) and write 

JV-l N 

« £ J>"G^(F fc £ - H - (b k £ - b k )D~U°M(x) 

3=1 k=l 
N 

= Y J D-G h ,{x,x k m - F° k - (b k - b k )D~U° k ). 
k=i 

Note that the expression above is an approximation because we have changed D~U £ on the 
right hand side of (3.8) to D~U°. The error is 

N 

rn(x) = -Y, D ~G%(x,x k )(b k £ - b k )D-{U e - U\. (3.22) 

k=l 

Estimating \D~ Gq| by Ch and using Cauchy-Schwarz on the sum over k and (3.13) and 
(3.12), we verify that E\r^{x)\ < Ce/r 5 /2||fl|| 1|H (i + ||/|| 2 ). 
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Using the expressions of 4> e and b £ , and the estimates of the higher order terms in them, 
(3.16), we can further approximate K\(x) by 



N 



X, Xk j 



k=l 



3>k-l 



t - Xk-l 

h 



q £ (s)ds 



Xk-l 



dt 



* 2 



a 



q e (t)dt 



%k-l 



The error in this approximation is: 
N / 



i 2 {x) = Y J D-G } ^ 



x,x k ) ri k 



k=i 



j-k 



fit) 



h 

%k-l 



Qe(s) 



q e (t)dt . 



t - X k -1 

h 



■i:k 



%k-l 



dt 



Using Lemma A.l, (3.17) and Cauchy-Schwarz, we have 



(3.23) 



E 



: (s)ds 



< Ce. 



Using this estimate, we verify that the mean of the absolute value of the first term in 
r\ 2 is bounded by Ce||/||2||-R||i,R. A similar estimate with \D~U®\ < Ch (due to super- 
convergence) shows that the second term in r\ 2 has absolute mean bounded by Cg:^ 1 1 1 .f? 1 1 i,m - 
Therefore, we have E|r^ 2 (a;)| < Ce/i _1 (l + ||/||2)||-R||i,r- We remark also that in the case of 
long range correlations, we should apply Lemma A. 2 instead. 

Moving on to K 2 (x), we observe that for fixed x, K 2 {x) reduces to a sum over at most 
two terms, due to the fact that 4>{ and (p J have local support only. Let j(x) be the index so 
that x G (xjf x )_i, Xjt x \]. We have 

N 

K 2 {x) =J2D-V}(4i - 4>l)(x) = D-C/J^M - ~4 x) ){x) 



q £ (t)dt 



x — x 



i(«0-i 



h 



j(x)-l 



: {t)dt 



In the second step above, we used the decomposition of 4> £ again. The error we make in 
this step is 



a*D-U° 
r h 2 (x) = r l3{x) ^ 



,(t)dt 



x — x 



j(x)-l 



h 



: (t)dt 



(3.24) 



x j(x)-l 



We verify again that E|r2(x)| < CeH-RHijR. 

Now for Kz(x), we use Cauchy-Schwarz and have 



N-l 



E\K 3 (x)\<E[ Y,\ U i- U l 



2 (N-l 



.3=1 



(3.25) 
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The last inequality is due to (3.12) and (3.10). 

In the approximations of K\{x) and K2(x), we change the order of summation and 
integration. We find that K\(x) is then f L\{x,t)q £ (t)dt plus the error term + r\ 2 , 
and K2{x) is Jq L2(x,t)q £ (t)dt plus the error term r^. Therefore we proved (3.18) with 
r £ (x) = r\i + r\ 2 + r \ + Ks(x). The estimates above for these error terms are uniform in 
x, verifying (3.19). □ 

3.3 Weak convergence of the corrector of a multiscale scheme 

In this section, we characterize the limit of the corrector — Uq, with proper scaling, in 
the multiscale scheme when e is sent to zero. As we have seen, the scaling depends on the 
correlation range of the random media. 

Proposition 3.3. Letu £ be the solution to (1.1) given by a multi-scale scheme that satisfies 
(N1)-(N3). Suppose (3.9) holds. Let be the standard FEM solution to (1.2). 

(i) Suppose that q(x) satisfies (S1)-(S3) and that the conditions of item one in Proposi- 
tion 3.2 hold. Then, 



V £ e_K) Jo 

(ii) Suppose that q(x) satisfies (L1)-(L3) and that the conditions of item two in Propo- 
sition 3.2 hold. Then, 



The real number a is defined in (2.3) and o~h is defined in Theorem 2.2. 

These results allow us to prove the weak convergence in step (ii) of the diagram in (1.4) 
for fairly general schemes. A standard method to attain such weak convergence results is 
to use the following Proposition; see [19, p. 64]. 

Proposition 3.4. Suppose {M £ } with e G (0, 1) is a family of random processes with values 
in the space of continuous functions C and M e (0) = 0. Then M £ converges in distribution 
to Mq as e — > if the following holds: 

(i) (Finite-dimensional distributions) for any < X\ < • • • < < 1, the joint distribu- 
tion of (M £ (xi), • • • , M £ (xk)) converges to that of (Mq(xi), • • • , Mo(xfc)) as e — > 0. 

(ii) (Tightness) The family {MeI^o,!) is a tight sequence of random processes in C{I). 
A sufficient condition is the Kolmogorov criterion: 35, (3, C > such that 



uniformly in e and t,s E (0, 1). 

The standard Kolmogorov criterion for tightness requires the existence of t G [0,1] and 
some exponent v so that sup e E|M e (f)| 1 ' < C for C independent of e and v. In our cases, 
since M £ (0) = for all e, this condition is always satisfied. 




(3.26) 




(3.27) 




(3.28) 
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We will prove item one of Proposition 3.3 in detail; proof of item two follows in the same 
way, so we only point out the necessary modifications. Recall the decomposition in (3.18). 
Let I £ denote the first member on the right hand side of this equation, i.e., the oscillatory 
integral. Let X h denote the right hand side of (3.26). The strategy in the case of SRC is to 
show that {e~2/ e } converges in distribution in C to the target process I h , while {s~^r £ } 
converges in distribution in C to the zero function. Since the zero process is deterministic, 
the convergence in fact holds in probability; see [5, p. 27]. Then (3.26) follows. 

Proof. Convergence of {e~?I £ }. We first check that the finite dimensional distributions of 
I £ (x) converge to those of I h (x). Using characteristic functions, this amounts to showing 

Eexp (± jf * q £ (t) E"=i $jL h (xj,t)dt \ ^ Eexp (ia jf * £?=i £ j L h &,t)dW t ) , 

for any positive integer n, and any n-tuples (x 1 , ■ ■ ■ ,x n ) and (£i, • • • , £ n )- We set m(t) = 
Yuj=i £jL \ xJ '' i • The convergence above is proved if we can show 

4= r qM)m(t)dt distribution ) a I' m(t)dW t , (3.29) 
V £ Jo 



o 



for any m(t) that is square integrable on [0,1]. Indeed, this convergence holds as long 
as q{x,w) is a stationary mean-zero process that admits an integrable p-mixing coefficient 
p(r) £ L 1 (R). This result is more or less standard, and a proof can be found in [2, Theorem 
3.7 and its proof]. Our assumptions (S1)-(S3) guarantee the existence of such a p(r). 
Therefore, we proved the convergence of the finite distributions of {e~2/ £ }. 

Next, we establish tightness of {e~^ I £ (x)} by verifying (3.28). Consider the fourth 
moments and recall L h = L\ + h\ in (3.20); we have 



E(I £ (x) - I £ (y)y < 8\E(^= I q e (t)(L^(x,t) - L^(y t t))dt 



q £ {t){L h 2 {x,t) - L h 2 (y,t))dt) 4 }. 



(3.30) 



We estimate the two terms on the right separately. For the first term we observe that 
Li(x,t) is Lipschitz continuous in x. This is due to the fact that G^x^Xk) is Lipschitz in 
x with a universal Lipschitz coefficient. Since the other terms in the expression of L\{x,t) 
in (3.20) are bounded by C, we have 

\L\{x,t) - L^{y,t)\ < -\x-y\. 
We use this fact and apply Lemma A. 4 to deduce 

E(-^= J q E (t)(Ll(x,t) - Ll(y,t))dty < ^\x - y\ 4 . (3.31) 
The constant C above depends on A, A, and ||p2|| 1K 
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To estimate the second term in (3.30), consider two distinct points y < x. Let j and k 
be the indices such that x G (xj-±,Xj] and y G (xk-i,Xk\- Then one of the following holds: 
j — k > 2, j — k = or j — k = 1. In the first case, since |-D~[/°| < Ch for some C depending 
on A, A and ||/||2, we have the following crude bound. 

\L h 2 {x,t)-L h 2 {y,t)\<C <^\x-y\. 

The same analysis leading to (3.31) applies, and the second term in (3.30) is bounded by 
C\x — y\ 4 /h 4 in this case. 

When \ j — k\ = 0, x and y are in the same interval (xj,Xj+\). We can write 



1 q £ (t)(L h 2 (x,t) - L h 2 (y,t))dt = — h ^ 

h 



q e {t)dt 



x-y 
h 



q £ (t)dt . (3.32) 



Since x and y are in the same interval, the function (x — y)/h is bounded by one. Now 
Lemma A. 4 applies and we see that the fourth moments of the members in (3.32) are 
bounded by 



C 



4f j\e{t)dt 



x-y 
h 



E( q £ (t)dt 



< C\x-y\ 2 . 



When j — k = 1, we have 



q £ (t)L^(y,t)dt 





( t- y 






a*D-U^ 


( t* 




\ Jy 



,{t)dt 



y-x~_ 2 f^- 1 



q e (t)dt 



h 

y - Xj-i 
h 



r(t)dt 



Xj-2 



q £ (t)dt 



Xj-2 



Let Xj-\ play the role of x in (3.32) and notice that L 2 (xj_i,t) = 0. We get 



1 



r-l \ 4 

-j= ^ q e (t)L^y,t)dt\ <C- 



y - Xj-i\ 



Similarly, in the interval where x lands, let play the role of y in (3.32). We have 



E(-L j\ £ (t)L h 2 (x,t)dtY < C l -^f 



We combine these estimates and see that in this case, the second term in (3.30) is bounded 
by 



l 



l 



8(E(^— q E (t)L%(y,t)dt) +E(^= J q E (t)L%(x,t)dt 
<c \ x i-i - y? + \x - xj~i\ 2 <c \x- y\ 2 



h 2 - h 2 ' 

In the last inequality, we used the fact that a 2 + b 2 < (a + b) 2 for two non-negative numbers 
a and b. 
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Combine these three cases to conclude that for any x,y £ [0,1], the second term in 
(3.30) is bounded by C\x — y\ 2 /h 2 . This, together with (3.31), shows 



E(-L^ q£ (t)(L h (x,t)-L\ y> t))dtY <C l -^f-. (3.33) 

In other words, {e~2/ e (x)} satisfies (3.28) with (3 = 4 and 5 = 1, and is therefore a tight 
sequence. Consequently, it converges to X h in distribution in C. 

Convergence of {e~2r^}. For the convergence of finite dimensional distributions, we 
need to show 

Eexp (V^X^iCM^)) 

for any fixed n, {x J }™ =1 and {£j}" =1 . Since \e %e — 1| < \6\ for any real number 9, the left 
hand side of the equation above can be bounded by 

-U Ei&^fo) <E;I&|4= sup E|r£(^)|. 

V £ V e l<i<" 

The last sum above converges to zero thanks to (3.19), completing the proof of convergence 
of finite dimensional distributions. 

For tightness, we recall that r^(x) consists of in (3.22), r^ 2 in (3.23), K^{x) in (3.18) 
and ri} (x) in (3.24). In the first three functions, x appears in Lipschitz continuous terms, 
e.g., in D~GQ{x;xk) or 4>i(x) — 4>q(x). Meanwhile, the terms that are x-independent have 
mean square of order 0(e) or less. Therefore, we can choose (3 = 2 and 5 = 1 in (3.28). For 
instance, we consider r^ 2 (x) in (3.23). Since q £ is uniformly bounded, the integrals of q £ on 
the interval Ik are bounded by Ch. Recall that \D~U%\ < Ch also; we have 

i 2 {x)-r^y)\<y^\D-{G h ^x,x k )-G h G (y,x k ))\(\f lk \ [ C\f\dt + C\r 2k 

By the Lipschitz continuity of D~Gq and the estimate (3.17), we have 

E ( rU*) -rUv) \ 2 < c l ]x _ y? supE{ | flfc | 2 + | f2fe | 2} < c \^yl_ 

\ V £ J £ k h 

Similarly, we can control r\-y and K3. For r\ in (3.24), we observe that it has the form 
of the main part of K2(x), which corresponds to L^x^t) and the second term in (3.30), 
except the extra integral of q e . Therefore, the tightness argument for the second term in 
(3.30) can be repeated. The extra q £ term is favorable: we can choose (3 = 2 and 5 = 1 in 
(3.28). 

To summarize, {e~2r^/y/e} can be shown to be tight by choosing (3 = 2 and 5 = 1 
in (3.28). Therefore, it converges to the zero function in distribution in C. We have thus 
established the convergence in (3.26). 

The case of LRC. In this case, the scaling is e™ . The proof is almost the same as above 
and we only point out the key modifications. 
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Let us denote the right hand side of (3.27) by X H . To show the convergence of the 
finite dimensional distributions of {e~^I £ }, instead of using (3.29), we need the following 
analogue for random media with LRC: 

i f 1 q e {t)m{t)dt distributi ° n > a n C m(t)dW t H , (3.34) 
£2 J e^O J 

where an is defined below (2.19). The above holds only for q(x,to) constructed as in (Ll)- 
(L3), and m £ L 1 (M) n L°°(R). This convergence result was established in Theorem 3.1 
of [3]. Assumptions on q(x) in item (ii) allows us to use this result and conclude that the 
finite dimensional distributions of {£~ ~ I £ (x)} converge to those of X H . For the tightness 
of {e~2l £ (x)}, we can follow the same procedures that lead to (3.31) and (3.32). We only 
need to consider second order moments when applying the Kolmogorov criterion thanks to 
Lemma A. 2, which says 

E (4r f V Qe(t)dty < C\x - y\ 2 ~ a . (3.35) 

In the SRC case, since a equals one we only have \x — y\ on the right. To get an extra 
exponent 5, we had to consider fourth moments. In the LRC case, a is less than one, so 
we gain a 5 = 1 — a from the above estimate. With this in mind, we can simplify the proof 
we did for (3.26) to prove that {£~2"/ e } converges to I H . Similarly, {e~^r £ } converges to 
the zero function in distribution, and hence in probability, in the space C. The conclusion 
is that (3.27) holds. This completes the proof of Proposition 3.3. □ 

From the proofs of the propositions in this section, the results often hold if the conditions 
in item (i) or (ii) of Proposition 3.2 are violated in an £-independent manner. For instance, 
if the second equation in (3.16) is modified to 

b k e - b k = c(h)[l + f 2k ] q e (t)dt) , (3.36) 

for some function c(h) and for region C Ik, then this modification will be carried to 
L h (x,t) and following estimates, but the weak convergences in Proposition 3.3 still hold. 



4 Weak convergence as h goes to 

In the previous section, we established weak convergence of the corrector u^ — Uq of a general 
multi-scale scheme when the correlation length e of the random medium goes to zero while 
the discretization h is fixed. In this section, we send h to zero, and characterize the limiting 
process. We aim to prove the following statement. 

Proposition 4.1. Let L h (x,t) be defined as in (3.20). As h goes to zero, the Gaussian 
processes on the right hand sides of (3.26) and (3.27) have the following limits in distribution 
in C: 

a C L h (x, t)dW t distributi ° n > U(x; W), (4.1) 
Jo 
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where IA is the Gaussian process in (2.14). Similarly, 

a H f 1 L h (x, t)dW t H distribution > U H (x; W H ), (4.2) 
Jo h -+° 

where Uh is the Gaussian process in (2.20). 

We consider the case of SRC first. Recall that I h (x) denotes the left hand side of (4.1). 
It can be split further into three terms as follows. Let us first split L\{x,t) into two pieces: 



MiOMj - 2^ lj fcW r ji — i 

k=l 

Lf 3 (x,t) = X;i/ k (t)a*D-G&(x,a: fc ) - / ' f(s)ds - - ' 

fc=l \ ^* •'^fc-l , 



(4.3) 



Then define Z^(x) by 

T\{x;W) = a I L h u {x,t)dW u i = 1,2. l£(x;W) = a [' L%(x,t)dW t . (4.4) 
Jo Jo 

As it turns out, X\(x\W) converges to the desired limit, while T^x^W) and X3 (x; W) 
converge to zero in probability. 

Proof of (4.1). Convergence of {T±(x)}. By Proposition 3.4, we show the convergence of 
finite distributions of {Zj^x)} and tightness. Since all processes involved are Gaussian, for 
finite dimensional distribution it suffices to consider the covariance function Ri(x,y) := 
E{X{ l (x)X{ l (y)}. By the ltd isometry of Wiener integrals, we have 

R h l {x,y)=a 2 [ L h n (x,t)L h u (y,t)dt. 
Jo 

For any fixed x, L^x, t), as a function of t, is a piecewise constant approximation of L(x, t). 
This is obvious from the expression of L(x,t) in (2.16). Therefore, L$i(x,t) converges to 
L(x,t) in (2.15) pointwise in t. Meanwhile, L\ x is uniformly bounded as well. The dominant 
convergence theorem yields that for any x and y, 

lim R^(x,y) = a 2 [ L(x,t)L(y,t)dt = E(U(x;W)U(y;W)). (4.5) 

h->0 Jo 

This proves convergence of finite dimensional distributions. 

The heart of the matter is to show that {Ui(x; W)} is a tight sequence. To this end, we 
consider its fourth moment 

4 r 4 4 
E(lf(x)-Xf(j/)) =/ \{{LUx,U)-LUy,tmX[dW u . (4.6) 



i=l i=l 
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Since increments in a Brownian motion are independent Gaussian random variables, we 
have 

4 4 

E [J dW u = Wi - ^(h - U) + 8{h - t 3 )5(t 2 - U) + 6(h - t 4 )5(t 2 - t 3 )] JJ dt { . (4.7) 

i=l i=X 

Using this decomposition, and the fact that the L\ x is piecewise constant, we rewrite the 
fourth moment above as three times of 

2 



1 \ 2 

(L h n (x,t) - L h n (y,t)) 2 dt 



o 



f a*D-(G%(x,x k ) - G%(y,x k ))a*D-U° k \ 2 

ItS h h ) 

Hence, we need to control HL^x, •)— £n(y, - ) 1 1 2 - Since Gq is the Green's function associated 
to (1.2), it admits expression (2.17) as remarked below Proposition 3.2. Fix y < x, and let 
ji and j2 be the indices so that y E (sc^-i, Xj-±] and x E (xj 2 -i, Xj 2 ]. Then we can split the 
above sum into three parts. In the first part, k runs from one to j\ — 1. In that case, both 
x k and x k -i are less than y. Formula (2.17) says: o,*(Gq(x, x k ) — G^(y,x k )) = x k (y — x). 
Consequently, 

a*D-(G%(x,x k )-G%(y,x k )) 



h 

since \D~U%/h\ is bounded, we have 



(y - x). (4.f 



<£' /«,-D-(Gt(x,xO-GS( y ,x t ))N' (^) 2 ft < C |x-/X> ^ CT-- „P. (4.9) 
fc=i ^ / V / fc=1 

Another part is A; running from j2 + 1 to A r . In that case, both x k and x k ~\ are larger than 
x. The above analysis yields the same bound for this partial sum. 

The remaining part is when k runs from j\ to j2- In this case, for some k, x k may end 
up in (y,x), and we have to use different branches of (2.17) when evaluating G^{x^x k ) and 
G(j(y, x k ). Consequently, the cancellation of h in (4.8) will not happen, and we need to mod- 
ify our analysis. We observe that, due to the Lipschitz continuity of Gq and boundedness 
of \D~U°/h\, we always have 



g f a*D~(G h (x,x k )-G h (y,x k )) \ 2 f a*D~U° k \ \ ^ < c \x - y\ 2 A ^ ^ 



If J2 — ii <• 1) the last sum above is then bounded by 2C\x — y\ 2 /h. In this case, it is clear 
that \x — y\ < 2h\ as a result, the sum above is bounded by C\x — y\. 

If j2 — ji > 2, the above estimate will not help much if j2~j\ is very large. Nevertheless, 
since \D~G\jh\ is bounded by some universal constant C . We have 

g f a*D-{G h {x,x k ) - G h (y,x k )) \ 2 / a'iTCff \ 2 ^ Q A ^ (4 n) 

fc=j'i ^ ^ fc=ii 
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Meanwhile, we observe that in this case 

3\x -y\> 3(x i2 _i - x h ) = 3(j 2 - ji - l)h = (j 2 - ji + l)h + 2(j 2 - ji - 2)h 
> {32-h + l)h. 

Consequently, the sum in (4.11) is again bounded by C\x — y\. Combining these estimates, 
we have 

\\L h n (x,-)-L h n (y,.)f 2 <C\x-y\. (4.12) 

It follows from the equation below (4.7) that {Z^(x)} is a tight sequence and hence converges 
to U(x,W). 

Convergence ofX^ to zero function. For the finite dimensional distributions, we consider 
the covariance function R 2 (x,y) = EjZ^x)^ (y)}- By ltd isometry, 



a 2 



L h 12 { Xl t)L h l2 (y,t)dt. (4.13) 



Now from the expression of L^ 2 (x, t), (4.3), we see that L^ 2 (x,t) converges to zero point- 
wise in t for any fixed x. Indeed, in the above expression, \D~G^/h\ is uniformly bounded 
while the integrals of f(s) and of /(s)^q(s) go to zero due to shrinking integration regions. 
Meanwhile, L\ 2 is also uniformly bounded. The dominated convergence theorem shows 
R 2 (x,y) — > for any x and y, proving the convergence of finite dimensional distributions. 
The tightness of {I 2 (x)} is exactly the same as {X^(x)}; that is to say, the properties of 
D~Gq can still be applied. We conclude that {X 2 (x)} converges to zero. 

Convergence ofl^(x) to zero. For the finite dimensional distributions, we observe that 
L 2 (x,t) is uniformly bounded and for any fixed x, it converges to zero point- wise in t, due 
to shrinking of the non-zero interval Ij( x y The covariance function of I^(x), therefore, 
converges to zero, proving convergence of finite dimensional distributions. 

For tightness, we consider the fourth moment of 1^{x) — Z^(y). By (4.7), it equals 
E(Z*(x; W) -ll(y,W)) A = 3 (^j\L h 2 {x,t)-L h 2 {y,t)fdt 



(4.14) 

Recalling the expression of L%(x,t) in (3.20), it is non-zero only on an interval of size h 
and is uniformly bounded. Let j(x) be the interval where L 2 (x) is non-zero, and similarly 
define j(y). Assume y < x without loss of generality. Consider three cases: j(x) = j(y), 
j(y) = j( x ) ~ 1; an d j(x) — j(y) > 2. In the first case, x and y fall in the same interval 
[xj_i,xj] for some index j. Then we have 

j\L h 2 (x,t) - L h 2 {y,t)fdt < C J 1 (l [Xjy] (t) - ^-^.(tjfdt. 
This integral can be calculated explicitly; it equals: 

"1 „ „, 



--(x -y)- 2^(x -y) + = (x - y)[l 



24 



Since |1 — (x — y)/h\ < 1 and \D U®/h\ < C, the above quantity is bounded by C\x — y\. 

In the second case, with j the unique index so that y < Xj < x and using the triangle 
inequality, we have 



\L h 2 {x,t) - L h 2 {y,t)\\l < 2 (\\L%(x,t) - L h 2 { Xj ,t)\\l + W^x^t) - L h 2 {y,t)\\ 



For the first term of the right hand side above, let Xj play the role of y in the previous 
calculation. This term is bounded by C(x — Xj). Similarly, for the second term, let Xj 
play the role of x, and we bound this term by C(xj — y). Consequently, we can still bound 
\\L*(x,-)-LZ(y,t)\\l by C\x-y\. 

In the third case, we have h < \x — y\. Meanwhile, since L 2 is uniformly bounded and 
is nonzero only on intervals of size h. We have 

\\L$(x,t)-L$(y,t)\\l<Ch<C\x-y\. 
Combining these three cases, the conclusion is: 

E(2*(x; W) - l^y, W)f < C\x - y\ 2 . (4.15) 
This proves tightness and completes proof of the first item of Proposition 4.1. □ 

In the proof above, we used the fact that Gq(x) defined in (3.21) is in fact the real 
Green's function defined in (2.17). However, the analysis follows as long as \D^Gq(x, Xk)/h\ 
is piecewise Lipschitz in x with constant independent of h, and the total number of pieces 
does not depend on h. 

The fact that 1 2 {x) and I3 (x) do not contribute to the limit is quite remarkable. It says 
the following. As long as the limiting distribution of the corrector — Uq is considered, 
the role of the multi-scale basis functions is mainly to construct the stiffness matrix, which 
is reflected by X^(x); its roles in constructing the load vector F 6 and in assembling the 
global function, which are reflected in T 2 (x) and 2g (x) respectively, are asymptotically not 
important. 

Now, we prove the second part of Proposition 4.1. 

Proof of (4.2). Recall that X H {x) denotes the left hand side of (4.2). Using the same 
splitting of L\ in (4.3), we can split T H into three pieces I Hi (x), i = 1, 2, 3, as in (4.4). The 
only necessary modification is to replace a by an and to replace the Brownian motion Wt 
by the fractional Brownian motion W t . We show that I H1 (x) converges to Uh while X^ 2 
and converge to the zero function. 

Convergence of finite dimensional distributions. For I H1 , we consider the covariance 
matrix R Hl (x,y) defined by M{X HX (x)X H1 (y)}. Using the isometry (A. 19), we have 

y) = k C t dtds . (4.I6) 



Jo 



As before, the integrand in the above integral converges to L(x,t)L(y, s)/\t 
every (t, s). Meanwhile, since L\ x is uniformly bounded, the integrand above is bounded by 
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C\t — s\~ a which is integrable. The dominated convergence theorem then implies that R H i 
converges to the covariance function oHAh{x] W h ). The convergence of finite distributions 
of I H2 and I H3 are similarly proved. 

Tightness. Due to LRC, we only need to consider the second moments in 3.28. For 
{Iflj}, we consider 

■eiiW -zjhfe)) 2 = » f (^"■ t )-^(»; t »(^(»-')-^(»-')) <M ., 

JR2 \t-s\ a 

using again the isometry (A. 19). Now we claim that 

\\L h ll (x,t)-L h ll (y,t)\\ LPt <C\x-y\K (4.17) 

for any p > 1. Indeed, for p = 2, this is shown in (4.12); the analysis there actually shows 
also that the above holds for p = 1. For p = oo, this follows from the uniform bound on 
L^. For other p, this follows from interpolation; see [21, p. 75]. 

Now, we apply the Hardy-Littlewood-Sobolev lemma [21, §4.3] to the expression of the 
second moment above. We obtain the bound 

C(a) K \\L^(x, •) - L h u (y, -)\\ Ll \\L h u (x, •) - L h n (y, .)|| i < C\x - y\ 2 ~ a . 

Therefore, the Kolmogorov criterion (3.28) holds with (3 = 2 and 5 = 1—a, proving tightness 
of {I hi}- Tightness of {^^2} foU° ws i n the same way because L\ 2 has the same structure 
as L\y as remarked before. Tightness of {2^ 3 } follows from the same argument above and 
the control on ||-L2( X > ') — ^(y, -)||| in the equation above (4.15). This complete the proof 
of (4.2). □ 



5 Applications to MsFEM in random media 

In this section, we prove Theorems 2.1 and 2.2 as an application of the general results 
obtained in the preceding two sections by verifying that the multiscale finite element method 
(MsFEM) satisfies the conditions of Proposition 3.3. 

5.1 The multiscale basis functions 

We describe the MsFEM for (1.1) as a special case of the methods developed in [16, 17]. 
We verify that this scheme satisfies assumptions (N1)-(N3). 

Recall that we have a uniform partition of the interval [0, 1] with nodal points {xk}^ =1 , 
where xo = 0, and xjy = 1, and 1^ the kth interval (xk-i, Xk)- The standard hat basis 
functions are denoted by {^(x)}^^ 1 . They span the space Vfr. The idea of MsFEM is to 

replace the hat basis functions by {4>i} which are constructed as 

C e 4(x) = o, xehul 2 u ...uitf_i, ^ 

= <Po> x& { x k}k=Q- 
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Here C £ is the differential operator in (1.1). Clearly, <fp e has the same support as 0q and 
thus satisfies (Nl). Note that the {</4} are constructed locally on independent intervals, 
and are suitable for parallel computing. 

For any k = 1, • • • ,N, we observe that the only non-zero basis functions are <p £ and 
0g -1 . Further, they sum up to one at the boundary points x k -\ and x k . Since equation 
(5.1) is of linear divergence form, we conclude that 4> £ (x) + 4> £ ~ l (x) = 1 on the interval. 
This shows that MsFEM satisfies (N3). In fact, the functions {4>e}k=i for MsFEM are 
constructed by solving (5.1) on I k with boundary values zero at Xk-i and one at x k - Once 
they are constructed, is given by (3.4). We can solve 4> £ analytically and obtain that 



V e T a~\t)dt, V e =( f a-\t)dt\ . (5.2) 



Consequently, (Nl) and (N3) indicate that MsFEM also satisfies (N2). To calculate the 
entries of the stiffness matrix we fix any i = 2, • • • ,N — 2, and compute 



(A^. u = - I a £ ^ \ dx = - ( a £ -^ ] / a^(s)ds = -b\. 




The last equality can be verified from the fact that 4>\ solves (5.1) and integration by parts. 
For i = l and N, we verify that (3.3) holds for b® and b^ given by (5.2). 

We record here a well-known super-convergence result: When dimension d = 1, the 
standard finite element method is super-convergent, in the sense that it yields exact values 
at nodal points. In our case, u^Xk) = uo(xk), where uq solves (1.2) and Uq is the FEM 
approximation. We observe that this property is preserved by MsFEM. Indeed, let Pu e 
be the projection of u £ in Vj 1 , i.e., Pu £ = u £ (xj)<pi(x). Then, using integrations by parts, 
(5.1), and the fact that Pu £ — u £ vanishes at nodal points, we have 



A £ (Pu £ , v) = A £ (u £ , v) = F(v), Vt; E V £ 



h 



Since the second equality is also satisfied by u £ , it follows that A £ (Pu £ — u £ ,v) = for any 
v in V £ . In particular, by choosing v = Pu £ — u £ and by coersivity of A e (-, •), we conclude 
that Pu £ = Ug. The super-convergence result follows. 

Several useful results follow from this super-convergent property. First, (x) of MsFEM 
coincides with the true solution u £ (x) at nodal points. Note that u £ can be explicitly solved 
analytically and that \u £ (x) — u £ (y)\ < C\x — y\ for some universal C. We then have 

\D-U%\ = \u h £ {x k ) - u h £ {x k ^)\ < Ch. (5.3) 

This improves the condition (3.9) in Proposition 3.2 and hence improves several subsequent 
estimates. Second, a fact which we have used extensively before, we have |-D _ Gq| < Ch and 
for any fixed x k , Gq (x; x k ) defined in (3.21) equals the continuous Green's function Gq(x, x k ) 
for (1.2). This is because the functions agree at the nodal points due to super-convergence 
and they are both piece- wise linear in x. 
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5.2 Proof of Theorems 2.1 and 2.2 

Since MsFEM is a scheme that satisfies (N1)-(N3), in order to apply (3.12) and (3.20) in 
previous propositions, we only need to check that (3.16) and (3.17) hold. 

Lemma 5.1. Let (f> £ and b k be the functions in (N1)-(N3) for MsFEM defined in (5.2). Let 
4>q and 6q be the corresponding functions for FEM. 

(i) Suppose a(x,uj) and q{x,u) satisfy (S1)-(S3). Then (3.16) and (3.17) hold and the 
conclusion of item one in Proposition 3.2 follows. 

(ii) Suppose a(x,uj) and q(x,uj) satisfy (L1)-(L3). Then the conditions and hence the 
conclusions of the second item of Proposition 3.2 hold. 

Proof. Prom the explicit formulas (5.2), we have 

-dt) -(! Ictt) =-b k £[ qe(t)dt. 



a e J \h k a* J h J±k 

Comparing with the second equation in (3.16), we find that it is satisfied with 

T2k ■= b k £ [ q £ (t)dt. (5.4) 



Similarly, we have 

r x k 




%(x) - = b k [ I q £ (s)ds - X X h k ~ l I q e (s)ds ) . (5.5) 



This shows again that (3.16) holds with fi£ having the same expression as fi\. defined 
above. In (5.4), since < b k < Ah~ l , we can apply Lemma A.l in the case of SRC or apply 
Lemma A. 2 in the case of LRC to conclude that E|f2fc| 2 < Ch~ 1 s in the first setting, while 
^l^2fc| 2 — C(eh~ l ) a in the second setting. This completes the proof. □ 

Note that the estimates (3.10) and (3.11) follow directly from this lemma. Therefore, we 
can apply Proposition 3.1 directly. Now we prove Theorem 2.1. Estimates (2.11) and (2.12) 
do not follow from Propositions 3.3 and 4.1 directly and need additional considerations. 

Proof of Theorem 2.1. Finite element analysis. We have seen that u £ super-converges to u £ . 
From (1.1) and (5.1), we observe that the following equation holds on Ij for j = 1, • • • ,N: 



C £ (u £ -u £ ) = f, in I. 
u £ — u £ = 0, on dij 



' J' (5-6) 



Using the ellipticity of the coefficient and integrations by parts, we obtain 



X\u £ - u e \ 2 H ij. < j a £ ^-(u £ - u £ ) ■ -^(u £ - u £ ) dx = J 



{u h £ - u £ )C £ {u h £ - u £ ) dx 



f{x){u h £ -u £ ){x) dx< \\f\\ 2 ,iM ~ u zhlr 
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Now recall that the Poincare-Friedrichs inequality says that 



l^-^lk/, < t-jII/IIv,, (5.8) 



he ~ U £ h,I 3 < -\u h £ - U £ \ H i j (5.7) 

Combining the inequalities above, we obtain 

\^-u £ \ mjj <±-\\f\\ 2Jj . 

Taking the sum over j, we obtain the first inequality in (2.11). To get the second inequality, 
we first apply the Poincare-Friedrichs inequality to the equation above to get 

Att 2 

and then sum over j. This completes the proof of (2.11) in item one of the theorem. 

Energy norm of the corrector. By energy norm, we mean the L 2 (fi, L 2 (D)) norm. Recall 
the decomposition of the corrector into Ki{x) in (3.15). For K\(x), we apply Cauchy- 
Schwarz to get the following bound for \K±\ 2 

EM - c/ ) 2 EiWto) 2 < EM - uf? (£i Ktxjf = \\u^- u% 2 . 

In the above derivation, we used the fact that (f^ (x) is non- negative, and ^ • <^q(x) = 1- Now 
we apply (3.12) to control this term. The function K 2 (x), as in the proof of Proposition 3.2, 
can be written as D~U^ x) x) - ^ 3 {x) ). Then from (3.10), we have E\K 2 (x)\ 2 < Ce||i?||i jK . 
For i^3, we have controlled E|i^3(x)| in (3.25). To control E\K 3 (x)\ 2 , we observe that 
1^3(^)1 < C || /|| 2- Note that all three estimates concluded in the three steps are uniform in 
x. Combining them, we complete the proof of (2.12). 

Convergence in distribution as e to zero. To prove item two of the theorem, we apply 
(3.26) of Proposition 3.3. We need to verify (3.9) in addition to (3.16) and (3.17), which we 
already verified in the previous lemma. But this is implied by (5.3), and hence we obtain 
(2.13). 

Convergence in distribution as h to zero. To prove (2.14), we apply the first result in 
Proposition 4.1. This completes the proof of the theorem. □ 

Proof of Theorem 2.2. In this case, the random processes q{x) and a{x) are constructed by 
(L1)-(L3). To prove the estimate in the energy norm, we follow the same steps as in the 
proof above, but use item two of Proposition 3.1 to control the term \\U e — U®\\ 2 2 in K\(x) 
and use Lemma A. 2 to control the terms in K 2 {x) and K%(x). 

To obtain the results in (2.19) and (2.20), we verify the conditions in item two of 
Propositions 3.3 and 4.1, applying the second case in Lemma 5.1 and following the steps in 
the previous proof. This completes the proof of the theorem. □ 



6 Applications to HMM in random media 

In this section, we adapt the general approach described in Sections 3 and 4 to the case of 
heterogeneous multiscale method (HMM). 
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6.1 The heterogeneous multiscale method 



The goal of HMM is to approximate the large-scale properties of the solution to (1.1) 
without computing the homogenized coefficient first. Suppose we already know this effective 
coefficient, i.e., a* in our case. Then the large-scale solution uq can be solved by minimizing 
the functional 



1 . , , „, , 1 f 1 , fdu\ 



1 /J„.\2 ,i 



I[u] := -A Q (u,u) - F(u) = - J a* I — j dx - j fu dx. 



o 



In numerical methods, the first integral above can be computed by the following mid-point 
quadrature rule: 



N / d \ 2 



Here x J = (xj-i + Xj)/2 is the mid-point of Ij. In HMM, a* is unknown, and the idea is to 
approximate (u'a*u')(x J ) by averaging in a local patch around the point x J . For instance, 
we can take 

■ du . -A Iff d(J>?u) 



Here, denotes the interval x- 7 + |(— 1, 1), that is, the small interval centered in Ij with 
length 5. The operator j£? maps a function to in V^, i.e., the space spanned by hat functions, 
to the solution of the following equation: 

(£ e (S?w) = 0, x e if u-.-u 

The idea here is to encode small-scale structures of the random media into the construction 
of the bilinear form. The key difference that distinguishes HMM and MsFEM is that the 
above equations are solved for HMM on patches it that are smaller than Ik- We check that 
«Sf is a linear operator; therefore, the following approximation of ^4o(') is indeed bilinear: 

N 



^Wj-Ej/,*^^*- (6.2) 

.7=1 I i 



3=1 j 

With this approximation of the bilinear form, HMM consists finding 

u^ ,S := argmin -A 5 £ (w,w) — F(w). 

This variational problem is equivalent to solving u^ ,S = U^' S (IPq(x), where U 6,5 is determined 
by the linear system 

A h,Sjje,5 = F 0_ 3) 

Therefore, the above HMM can be viewed as a finite element method. The finite dimensional 
space here is Vq . Therefore HMM clearly satisfies (Nl) and (N3). To check (N2), we 
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calculate the associated stiffness matrix As ,S . It has entries A £ ((f> z , <fii). From the defining 



equation (6.1), we see that is non-zero only on lfulf +l , which implies that A £ ' is again 
tri-diagonal. Further, we verify that J£4>q + =£?</>q -1 = 1 on the interval if, which can be 
obtained from integrations by parts and which implies that A £ ,S satisfies (3.3). Therefore, 
HMM satisfies (N2). 

In fact, we can calculate the b £ vectors. Let us consider the (i — l,i)th entry of A £ ' S , 
where i can be 2, • • • , N — 1. Since (J*f(/>Q -1 )' = — on if , we have 

A h A =- h -\ a e (s)(^l) 2 ds. 

) ei-li 5 JjS £ \ dx J 

Now from (6.1), we verify that a £ (^f(f> l )' on if is a constant given by 



C; 




Therefore, we have 

( Ak ' 5 ) £l _ u = ~^) 2 - 5 I . a e l ds = -~ ( / . a-Hs)ds ) =: (6.4) 

We extend the definition of b £ s to the cases of i = 1 and i = N, and check that the (1, l)th 
and (N — 1,N — l)th entries of the stiffness matrix also satisfy (3.3). In particular, the 
action of A £ ' S on a vector satisfies the conservative form as in (3.5). In the sequel, to 
simplify notation, we drop the 5 in the notations A £ ' , U 6,s and b % £ s . 

The well-posedness of the optimization problem above, or equivalently of the linear sys- 
tem (6.3) is obtained by Lax-Milgram. We show that the bilinear form Af (•, •) is continuous 
and coercive. Consider two arbitrary functions w = Wi4> l and v = VjCpQ in Vq. Then: 

A h E {w,v) = WiA^Vj = -ZiWiD+ibiD-Vi) = E^-Wi^IT^. 
Estimating the entries of vector b e by its infinity norm and use Cauchy-Schwarz, we obtain 

\A h £ {w,v)\ < ( sup bi) \\D^W\\ e 2\\D~V\\ e 2 < A\w\ H i \v\ H i . 

\l<i<N J 

In the last inequality above, we used the fact that A/i" 1 < b £ < Ah , which can be seen 
from its definition in (6.4) and the uniform ellipticity of a £ , and that ||-D - VF||^2 = \ w\jji\^h 
for w £ Vq. This proves continuity. Taking w = v , we have 

A h £ {w,w) > (w£ N b{\ \\D~W\\ e 2\\D-W\\ e 2 > X\w\ 2 Hl . 

This proves coercivity. Therefore, by the Lax-Milgram theorem for the finite element space 
[24, p. 137], there exists a unique Ue' S £ Vq that solves the optimization problem. Further, 
we have 

\u h /\ m <\ sup J-||/|| 2 . (6.5) 

An immediate consequence is that \D~~U e \ < C\fh by the argument following Proposition 
3.1. 



31 



6.2 Proof of Theorem 2.3 



To prove Theorem 2.3, we apply Proposition 3.2 to write the corrector u^' 5 — Uq as an 
oscillatory integral plus a lower order term. To apply Propositions 3.3 and 4.1 and obtain 
the weak convergences, we need to consider the difference b k — 6q since <f>i = in HMM. 
From the expression of b k £ in (6.4), we have 

b k £ ~b k = -b k e ^ I q e (t)dt = -(1 + fi»)~ f q E (t)dt, 

where f<i k is a random variable defined by 



r 2 k = -7&* / q e {t)dt 
6 H 



We verify that in the case of SRC, i.e., when q{x) satisfies (S1)-(S3), we have 



F f \ 2 F 

|2 «^ fi II oil w uk uk* 



n? 2k \ 2 < C-\\R\\ ltR , E(&£-6gJ <C^||i2||i, R , (6.6) 

for some universal constant C. Comparing this with (3.11) and (3.17), we observe that the 
estimates have been multiplied by a factor 4 in the HMM case. Similarly, it can be checked 
that in the case of LRC, i.e., when q{x) satisfies (L1)-(L2), these estimates will be multiplied 
by a factor of (4) ■ With these formulas at hand, we prove the third main theorem of the 
paper. 

Proof of Theorem 2.3. Short range media and amplification effect. In this case, the differ- 
ence of b £ — 6q and an estimate of it was captured in (6.6) and the equation above it. We 
cannot apply Propositions 3.2 and 3.3 directly. However, due to the remark at the end of 
section 3.3, similar conclusions still hold. The same procedure as in the proof of Proposition 
3.2 shows that the L h {x,t) function for HMM is: 

M-sl^hfi) 1 1 — • ( 6 - 7 ) 

k=l 

The first weak convergence in (2.21) holds with this definition of L h,s as an application of 
a modified version of Proposition 3.3. Indeed, the proof there works with L h ' S playing the 
role of L-q. The tightness is still obtained from the function D~Gq, and the factor h does 
not play any role at this stage. 

When h goes to zero, we can follow the proof of Proposition 4.1 to verify the second 
convergence in (2.21). Indeed, tightness can be proved in exactly the same way. All that 
needs to be modified is the limit of the covariance function of U h,s (x; W), which is defined 
to be a J 1 L h < 5 (x, t)dW t . This covariance function, by the Ito isometry, is as follows: 



R h ' 6 (x,y) :=a 2 C L h > 5 (x,t)L h ' S (y,t)dt 
Jo 

2 h 2 ^ a*D-G%(x;x k )a*D-G%(y-,x k ) fa*D~U% 

= a P^ 6 h h 

k=i 



(6.8) 
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Recall the expression of in (4.3). We verify that the above quantity can be written as 

a 2 - 8 j\ h n (x,t)L^(y,t)dt. 

Now the convergence in (4.5) implies that R h ' S converges to the covariance function of 
J^U(x; W). This completes the proof of (2.21). 

Long range media. The expression for — 6q are given above. Therefore, we can apply 
Propositions 3.2 and 3.3 (with modifications), to show that as e goes to zero while h is 
fixed, the HMM corrector indeed converges to (x;W H ) defined in (2.22). When h is 
sent to zero, we can follow the proof of Proposition 4.1 and show that U^f converges in 
distribution to some Gaussian process. To find its expression, we calculate the covariance 
function of . Thanks to the isometry (A. 19), it is given by 

<(*■»):=»/' f ^^h n, (6.9) 
Jo Jo l r — s \ 

Using the expression of L h ' S , and the following short-hand notations: 

j , , a*D-G%(x;x k )a*D-U° k 
J k{x ) := _ 

the covariance function can be written as 

nh 2 / rT/ , T/S T/NT/N1 /" f dtds . -v T / s -r r s f f dtds 



f f dtds 

[Jk{x)J m (y) + J m (x)J k (y)] / — + Jk( x )Jk(y) 

hs Its \t — s\ a 



o hs hs \t — s\ a hs Its \t 

k=lm=l 'k 1 ™ 1 1 fe=l k k 1 

The integral of \t — s\~ a can be evaluated explicitly: 

N fc-1 , 9 



£ Yl i J k(x)J m (y) + J m (x)Jk(y)]-^(i(k - m)h + 5} 2 ~ a + 



v ;v ; fc=lm=l 

N h 2 

-2p - m)h] 2 ~« + [(* - m)h - 5] 2 ~ a ) + , £ 2 J k (x)J k (y)-5 2 - a . 

J (l-a)(2-a)^ 5 2 

When m < k, the quantity between parentheses together with the 5 2 on the denominator 
forms a centered difference approximation of the second order derivative of the function 
r 2 ~ a , evaluated at {k — m)h, i.e., at t — s. This derivative is precisely (1 — a)(2 — a)\t — s\~ a . 
Meanwhile, the h 2 on the nominator can be viewed as the size of the measure dtds on each 
block I k x I m . Furthermore, Jk{x) is precisely L > l 1 (x,t) evaluated on I k . The conclusion is: 
those terms in the above equation with m < k form an approximation of 

o Jo \t — s \ a 



The second sum corresponds to the diagonal terms k = m. Since \J k \ is bounded, this sum 
is of order 0{h5~ a ) and does not contribute in the limit as h — > 0, as long as 5 3> ha . R h ^ 
converges to the covariance function of Uh(x\ W ), finishing the proof of (2.22). □ 
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7 A hybrid scheme that passes the corrector test 



We now present a method that eliminates the amplification effect of HMM with 8 < h 
exhibited in item one of Theorem 2.3 when the random media has SRC. Such an effect 
arises because the fluctuation in the short-range averaging effects occurring on the interval 
of size h are not properly captured by averaging occurring on an interval of size 5 < h. 

The main idea is to subdivide the element Ik uniformly into M smaller patches and 
perform M independent calculations on each of these patches. This is a hybrid method 
that captures the idea of performing calculations on small intervals of size 5 <C h to reduce 
cost as in HMM while preserving the averaging property of MsFEM by solving the elliptic 
equation on the whole domain. 

Let 6 = h/M be the size of the small patch i| for 1 < £ < M. Define b^ by 



(7.1) 



This definition is motivated by (6.4). Given a function w in the space V h , we define its 
local projection into the space of oscillatory functions in the small patches i| by: 




s€u£L 



i u e= 



f)T e 
1 aJ fc; 



(7.2) 



where w\ denotes this local projection. Recall that 4>q is the left piece of the hat basis 
function. Integrations by parts show that b^ = A e (((f>Q)e, (</>o)^), where A £ is the bilinear 



form defined in (3.2). HMM choose one small patch I,* and uses A e {(^)q)^, (0q) 
to approximate the value ^4 £ (0q,0q). Of course, the scaling h/5 is needed. This scaling 
factor turns out to amplify the variance as h goes to zero when the random medium has 
SRC. 

We modify the method of HMM by constructing b £ as b\ := ^2dLi b k eV In other words, we 



define the stiffness matrix element A^_ u by — YleLi ^e((0o 
we verify that 



With this definition, 



M 



a*- x ds 



1=1 L k J1 k 



{s)ds 



l ds 



Rewriting the sum of the first terms in the parenthesis, we obtain 



where r £ accounts for the sum over the second terms in the parenthesis. Clearly, E|r^| < 
Ce(M) -1 . This decomposition of b £ — bo and the estimate of r £ shows that we can apply 
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Proposition 3.2 to obtain the decomposition of the corrector. The L h (x,t) function in 
this case will be Lji(a;, i) in (4.3). Then it follows from Propositions 3.3 and 4.1 that the 
corrector in this method converges to the right limit. 

In this modified method, all the local informations on 1^ are used to construct b\ as in 
MsFEM. The main advantage is that the computation on {ifypL^ can be done in a parallel 
manner. The calculation in MsFEM performed on a whole domain of size h is replaced by 
h/S independent calculations. Accounting for the coupling between the h/5 subdomains is 
necessary in MsFEM. It is no longer necessary in the modified method, which significantly 
reduces is complexity. 

8 Conclusions and perspectives 

This paper analyzes the theory of random fluctuations for several multi-scale algorithms 
such as MsFEM and HMM in the simplified setting of a one-dimensional elliptic equa- 
tion. One of our main results is that MsFEM and HMM with 5 = h correctly capture the 
random fluctuations beyond the homogenization limit when the random media satisfy the 
short-range correlation (SRC) assumptions or the long-range correlation (LRC) assump- 
tions considered in this paper. Such schemes pass the corrector test we have introduced 
and are therefore likely to be reliable in the analysis and the assessment of uncertainty 
quantifications even in settings where a well-understood theory of random fluctuations is 
not available. 

We have also shown that less expensive schemes such as HMM with 5 < h still capture 
the random fluctuations in media with LRC but may fail to capture them in media with SRC. 
The reason is that media with LRC display self-similarities across scales. As a consequence, 
the properties of the random fluctuations are correctly captured at the macroscopic scale h. 
On the other hand, media with SRC display averaging effects at the scale e <C h. Numerical 
schemes thus need to solve the equation at the microscopic scale on the whole domain in 
order to capture such effects. This forces the choice 5 = h in the standard HMM scheme 
and led us to the modification of the MsFEM and HMM schemes presented in section 7 to 
pass the corrector test for both LRC and SRC media. 

The random fluctuations considered here all have a Gaussian structure and follow from 
an application of the central limit theorem for SRC media and a different but similar 
averaging procedure for LRC. Random fluctuations corresponding to rare events have to be 
analyzed with different methods. We refer the reader to [4] for the one dimensional elliptic 
equation in the setting of large deviation theory. 

The theory presented in this paper is restricted to the one-dimensional case. The main 
reason is that the theory of random fluctuations is well understood only in one dimension; 
see [3, 6]. In this setting, the random elliptic solution can be written as a stochastic integral 
whose analysis follows from fairly standard techniques. Note that the discretized solution 
solves a discrete system whose inversion is not explicit. This explains in large part why the 
results presented in this paper are rather long and technical. 

In higher dimensions, only specific equations are amenable to fluctuations analyses; see 
[2, 14]. For such equations, the homogenization limit is trivial. Unlike the case considered in 
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this paper, purely deterministic discretizations and basis functions based on the unperturbed 
elliptic equation provide accurate numerical solutions. Such models thus do not capture the 
difficulties inherent to the simulation of elliptic equations. They do not allow us to analyze, 
for instance, the resonance effects observed for MsFEM schemes and the over-sampling 
techniques designed to eliminate such effects [12, 16, 17]. The theory of fluctuations and its 
effect on multi-scale algorithms remains an open problem in such a setting. 

Even though the one dimensional equation is of purely academic interest, the analyses 
presented in this paper quantify the influence of using less expensive numerical schemes 
on the accuracy of the numerical solutions. In the case of SRC, random fluctuations are 
captured only by fairly expensive schemes that solve the elliptic equation on the whole 
domain of interest. In the case of LRC, self-similarity properties allow us to use much 
coarser schemes without sacrificing accuracy. These guiding rules should prevail for a large 
class of equations, elliptic and non-elliptic, linear and nonlinear, and hopefully in practical 
settings where a well-controlled theory of random fluctuations is not available. 
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A Several useful lemmas 
A.l Two moments 

We first record two key estimates on oscillatory integrals, which we have used already. The 
first one accounts for random media with SRC. 

Lemma A.l. Let q(x,uj) be a mean-zero stationary random process with integrable corre- 
lation function R(x) . Let [a, b] and [c, d] be two intervals on W and assume b — a < d — c. 



Proof. Let T denotes the expectation of the double integral. It has the following expression: 
T= ! ! R( t -^)dtds = [ [ )l M (t)l M ( S )dtda. 

J a Jc £ JRJR e 

We change variables by setting t —¥ t and (t — s)/e — > s. The Jacobian of this change of 
variables is e. Then we have 



Then 




(A.l) 



T\ < e 



JR JR 



[ [ \R(s)\l [aM (t)dtds = e(b - a)\\R\\ 1M 



This completes the proof. 



□ 
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The second one accounts for a special family of random media with LRC. The proof is 
adapted from [3]. 

Lemma A. 2. Let q(x,cj) be defined as in (L1)-(L3). Let F be a function in the space 
L°°(M). Let (a, b) and (c, d) be two open intervals and assume b — a < d — c. Then we have 



E 



b r d 



qC-)q{- £ )F{t)F{s)dtds 



< Ce a (b-a)(d-c) 



l-a 



(A.2) 



The constant C above depends only 071 K, (X CLTld Halloo, R* 
Proof. By the definition of the correlation function R, we have 



f j%{ t -)q{ S -)F(t)F(s)dtds} 



e-«R(l^)F(t)l [aM (t)F(s)l [c ^s)dtds. 

E 



As shown in [3], R(t) is asymptotically kt a with k defined in (2.8). We expect to replace 
R by KT~ a in the limit. Therefore, let us consider the difference 



e -«i2(L_i) _ 

v e 1 t-s a 



\F(t)\l [aM (t)\F(s)\l [c4] (s)dtds. 



By the asymptotic relation R ~ KT~ a , we have for any 5 > 0, the existence of T$ such that 
\R(t) — KT~ a \ < 5r~ a . Accordingly, we decompose the domain of integration into three 
subdomains: 

£>i = {(t, s ) e R 2 ,\t- s\ < T 5 e}, 
D 2 = {(t, s) £ M 2 , T s e <\t-s\< 1}, 
D 3 = {(t,s) £K 2 ,1< \t-s\}. 
On the first domain, we have 



e- a R 



t - s, 



K 



< 



D, 



e ' \t-s\ a 
,t- s 



\F(t)\l [aM (t)\F(s)\l [c4] (s)dtds 



e~ a R(- 



£ 



Fi(t)||F 2 (s)|dtda + 





K 


L 


\t - s\ a 



|Fi(t)||F 2 (s)|dfofo. 



Here and below, we use the short hand notation F\ (t) = F(t)l[ at b](t) andF 2 (s) = F(s)l[ c d i(s). 
The above integrals are then bounded by 



R\\oo,m [ \F(t)\ [ + SE \F 2 (s)\dsdt + / \F(t)\ f 5£ K\s\- a \F 2 (t - s)\dsdt 

J a Jt-T s e J a J-Tze 

2kT. 



<\\F\ 



S,oo 



1-ct. 



1 — a 



l-a 



On domain D 2 , we have 

,t — S, K 



D 2 



£~ a R{- 



£ 



t - S 



|Fi(t)||F 2 (s)|dtds < <5 / |t-s|^ a |Fi(t)||F 2 (s)|dtds 
2S\\F\ 



pb pi 2(5111 

<2S \F(t)\ \s\~ a \F 2 (t - s)\dsdt < ^2^ ( 1 + r i-V- Q ). 

J a JT s e 1 
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On domain D 3 , we can bound \t — s\ a by one, and we have 



D 3 



K £ ' t-S 



|Fi(t)||F 2 (s)|dtds < 5 / |t-sr Q |Fi(t)||F 2 (s)|(it(is 
Jd 3 



<2S [ b \F(t)\ f 1 \F(t-s)\dsdt<25\\F\\ 2 ^ R (l + T s e). 

Ja JT s e 

Therefore, for some constant C that does not depend on e or 5, we have 

»6 fd 



lim sup e 



E 



la J c 

Sending S to zero, we see that 

,-b r d 



(i) q ( S -) Fl (t)F 2 (s) 



t - S\ 



Fi(t)F 2 (s] 



< C\\F\\' R 6. 



lim e~ a E 



q^q^F^F^dtds 



-F l {t)F 2 {s)dtds. 



\t - s\ a 

Finally, from the Hardy-Littlewood-Sobolev inequality [21, §4.3], we have 



F 1 (t)F 2 (s) 



dtds 



\t - s\ a 
This completes the proof. 



< CUFxIli^ll^ll^-i^ < C||F||^(6 - a)(d - c 



l-Q 



(A.3) 

(A.4) 
□ 



A. 2 Moment bound for stochastic process 

In this section we provide a bound for the fourth order moment of q(x,co) in terms of the 
L 1 norm of the p-mixing coefficient. 

Let V be the set of all ways of choosing pairs of points in {1, 2, 3, 4}, i.e., 

V := [p= {{p(l),p(2)},{p(3),p(4)}} I KOe {1,2,3,4}}. (A.5) 
There are C| = 15 elements in V. 

Lemma A.3. Let q(x,u) be a stationary mean-zero stochastic process. Assume E|g(0)| 4 is 
finite and q(x,uj) is p-mixing with mixing coefficient p(r) that is decreasing in r. Then we 
have 

4 

|E{JJg(^)}| <E|g(0)| 4 ^pl(|^ (1) -x K2) |)pl(|x K3) -x K4) |). (A.6) 
i=i P ev 

Proof. Given four points {q(xi)}, i = 1, ■ ■ ■ , 4, we can draw six line segments joining them. 
Among these line segment there is one that has the shortest length. Rearranging the 
indices if necessary, we assume it is the one joining x\ and x 2 . Then set A = {xi,x 2 } and 
B = {x3,x^}. Rearranging the indices among each set if necessary, we assume also that 
d(A, B) is obtained by \x\ — x$\. Then by the definition of p-mixing, we have 

4 

|E{JJg(a; i )} - R(xi - x 2 )R(x 3 - x 4 )| < Var{q(xi)q(x 2 )} 2 Var{q(x 3 )q(x 4: )} ? p(\xi - x 3 \). 
i=l 

(A.7) 
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We can bound Vem:{q(xi)q(x 2 )} , and similarly the variance of \ai{q{x^)q{xi)} , from above 
by (E|g(a;i)| 4 E|g(x2)| 4 ) 1/2 . Therefore, the expression above can be bounded by E|(j(0)| 4 p(|:ri- 

*3|)- 

Since p is decreasing and \x\ — x%\ > \x\ — x 2 \, we also have 

4 

lE-tn^i)} - R ( x i - x 2)R{xz - x A )\ < E\q(0)\ 4 p(\ Xl - x 2 \). (A.8) 
i=i 

Now observe that min{a,6} < (ab)* for any two non-negative real numbers a and b. Ap- 
plying this observation to the bounds of the two inequalities above, and using the triangle 
inequality, we obtain 

4 

|E{U«(a=i)>| < \R(xi-x 2 )\ ■ \R(x 3 -x A )\ + E| g (0)|V(l*i ~ x 2 \)pH\ Xl "^D- (A.9) 

i=l 

Using the definition of mixing again, we obtain that \R(xi — x 2 )\ = \Eq(xi)q(x 2 )\ is bounded 
by 

yaA(q(x 1 ))\aA(q{x 2 ))p(\x l - x 2 \) < (E|g(0)| 4 )^ (\ Xl - x 2 \). 

In the last step, we used the fact that p < p? since p can always be chosen no larger than 
1. We can bound R{x\ — x^) in the same way. Therefore, we obtain 

4 

< %(0)| 4 [pH\xi - x 2 \)pi{\x 3 - x 4 \) + pH\ Xl - x 2 \)p*(\xi - x s\)} ■ (A.10) 

This completes the proof. □ 

We now derive a bound for the fourth order moment of oscillatory integrals of q e . 

Lemma A. 4. Let q{x,uj) satisfy the conditions in the previous lemma. Assume in addition 
that the mixing coefficient satisfies that ||/92||ir, is finite. Let (x,y) be an interval in R. 
Then for any bounded function m(t), we have 



E 



j= J\{^)m(t)dt\ <60E|g(0)| 4 -||^|| 2 iR JH|^-|x-y| 2 . (A.ll) 



Proof. The left hand side of the desired inequality is 

4 , 4 



i ry ry rv ry ' +. JL. 
/=-/ / / / EHqC^HmiuntthtsU}. (A.12) 

£ Jx Jx Jx JX ■ 1 ^ -_i 



i=l i=l 



Here and below, d\t\ ■ £4] is a short-hand notation for dt± ■ ■ ■ d A . Apply the preceding lemma, 
we have 



PGP' 



X J X Jx J X 
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Note that we did not write absolute sign for the argument in the p functions. We assume 
p is extended to be defined on the whole R by letting p{x) = p{\x\). There are 15 terms in 
the sum above that are estimated in the same manners. Let us look at one of them, with 
p{\) = p(3) = 1, p(2) = 2, and p{4) = 3. We perform the following change of variables: 

t\ — t2 t\ — t 3 

> t2, > ts, t± — > t±, £4 — > £4. 

e e 

The Jacobian resulting from this change of variable cancels e 2 on the denominator. The 
integral becomes 

t^—X t-^—x 



dti / dU / P*{t2)dt 2 / pz{t 3 )dt 3 . (A.13) 
J x Jb-x Jh-y. 

£ £ 

This integral is finite and is bounded from above by 

\x ~ y\ 2 \\pHl,n = 4|x-y| 2 ||p5|| 2 )R+ . (A.14) 
The other terms in the sum have the same bound. Hence we have, 

/ < E|g(0)| 4 x 15 x A\x - y\ 2 \\p^ || 2 , K+ ||m||^. (A.15) 
This verifies (A.ll) and completes the proof. □ 



A. 3 Fractional Brownian motion 



For the convenience of the reader, we briefly review some essential properties of fractional 
Brownian motion (fBm), and stochastic integral with fBm integrator. 

A fBm W H (t) with Hurst index H is a mean-zero Gaussian process with W H (0) = 0, 
stationary increments and ff-self-similarity, that is, for a > 0, 

{W H (at)} tm 1 {a H W H (t)} tm , (A.16) 

where = means the equality in the sense of finite dimensional distributions. From this simi- 
larity relation, we deduce E[(W H (t)) 2 } = \t\ 2H E[(W H (l)) 2 }. In particular, if E[(W H (1)) 2 } = 
1 we say the fBm is standard. It follows from the stationarity of increments that the co- 
variance function of W H (t) is given by 

R H (t, s) = E{W H (t)W H {s)} = i (\t\ 2H + \s\ 2H -\s- t\ 2H ) . (A.17) 

When H = 1/2, the increments of W H are independent and the fBm reduces to the usual 
Brownian motion. For H ^ 1/2, the increments are stationary but not independent. 

Stochastic integrals with respect to fBm can be defined on many functional spaces. Note 
that H = 1 — % is in the interval (^,1). In this case, a convenient functional space to define 
stochastic integral is 

\T\ H = {f-- f f |/(x)||/(y)||x- 2 /| 2 ^- 1 )^y<ool. (A.18) 
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It is easy to check, for instance from Hardy-Littlewood-Sobolev lemma [21, §4.3], that 
L 1 (M) n L 2 (R) C L x l u C |r|^. Stochastic integrals against fBm do not satisfy Ito isometry; 
instead, we have 

E { / f(t)dW t H [ h(s)dW s H ) = H(2H - 1) / u m 5^ H) dtds. (A.19) 
Urn Jr J Jm? \t - s] 2 ^-") 

For a nice review on stochastic integral with respect to fractional Brownian motion, we 
refer the reader to [23]. 
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